(load "infix")

(defmacro funcall. (f &rest args)
  `(the double-float
    (funcall ,f ,@(loop for arg in args
			collect `(the double-float ,arg)))))

(defun midpoint (f a b slices)
  "Integrate f (a function of one double-float which produces another
double-float) from a to b where a and b are both double-floats. Uses
the midpoint method with a given number of slices where slices is a
fixnum."
  (declare (double-float a b)
	   (fixnum slices)
	   (function f))
  (let* ((dx (/ (- b a)
		slices))
	 (current-point (+ a (/ dx 2))))
    (loop for n below slices
	  sum (* (the double-float
		   (funcall f current-point))
		 dx)
	  do (incf current-point dx))))

(defun simpson (f a b slices)
  "Integrate f (a function of one double-float which produces another
double-float) from a to b where a and b are both double-floats. Uses
Simpson's rule with a given number of slices where slices is a
fixnum."
  (declare (double-float a b)
	   (fixnum slices)
	   (function f))
  (let* ((dx (/ (- b a)
		slices))
	 (x a)
	 (sum 0d0))
    (loop for i from 1 to (1- slices) by 2
	  do (progn (incf sum (funcall. f x))
		    #I(x = a + i*dx)
		    (incf sum (* 4d0 (funcall. f x)))
		    (incf x dx)
		    (incf sum (funcall. f x))))
    #I(dx/3 * sum)))

(defun test-fun (x)
  (declare (double-float x))
  (the double-float
    #I(3*x^^2 + 2)))

;; int[3*x^2 + 2,x,2,5] = 123
(time (midpoint #'test-fun 2d0 5d0 500000))
(time (simpson #'test-fun 2d0 5d0 10000))


(defun euler (f x0 y0 &key (step 0.1d0)
	                   (steps 30)
	                   (stream *standard-output*))
  "Use Euler's method to output the curve of a differential equation given
   initial values of t and y(t) (think of them as x0 and y0)."
  (declare (function f))
  (loop for count below steps
	for x = x0 then (+ x step)
	for y = y0 then (+ y (* step
				(funcall. f x y)))
	do (format stream "~F  	~F~%"
		   (coerce x 'single-float)
		   (coerce y 'single-float))))

(load (translate-logical-pathname "clocc:src;cllib;simple"))
(load (translate-logical-pathname "clocc:src;cllib;gnuplot"))
(use-package :cllib)

(eval-when (:compile-toplevel :load-toplevel :execute)
  (setf cllib::*gnuplot-path* "/usr/bin/gnuplot")
  (setf cllib::*gnuplot-msg-stream* (make-broadcast-stream)))

(defun plot-euler (f x0 y0 &key (step 0.1d0) (steps 30))
  ;; FIXME: this temporary file generation is lame
  (let ((temp-filename (format nil "/tmp/euler-gnuplot-~A.tmp"
			       (get-universal-time))))
    (with-open-file (temp-file temp-filename
			       :direction :output
			       :if-exists :supersede)
      (with-plot-stream (gnuplot-stream :plot :plot)
      	(euler f x0 y0
	       :step step
	       :steps steps
	       :stream temp-file)
	(format gnuplot-stream "set terminal x11~%")
	(format gnuplot-stream "plot '~A' using 1:2 notitle with lines~%"
		temp-filename)))
    (delete-file temp-filename)))

;; dy/dt = 18y^2
(defun diffeq2 (x y)
  (declare (ignore x)
	   (double-float x y))
  (* 18d0 y y))

(plot-euler #'diffeq2 -2d0 -0.1d0
	    :step 0.01d0
	    :steps 3000)

(with-open-file (f "euler2.dat"
		   :direction :output
		   :if-exists :supersede)
  (euler #'diffeq2 -2d0 -0.1d0
	 :step 0.01d0
	 :steps 300
	 :stream f))

(defun diffeq (x y)
  (declare (double-float x y))
  (the double-float
    (- x y)))

(with-open-file (f "euler3.dat"
		   :direction :output
		   :if-exists :supersede)
  (euler #'diffeq -1d0 4d0
	 :step 0.01d0
	 :steps 300
	 :stream f))